R Notebook

The following is your first chunk to start with. Remember, you can add chunks using the menu above (Insert -> R) or using the keyboard shortcut Ctrl+Alt+I. A good practice is to use different code chunks to answer different questions. You can delete this comment if you like.

Other useful keyboard shortcuts include Alt- for the assignment operator, and Ctrl+Shift+M for the pipe operator. You can delete these reminders if you don’t want them in your report.

---

Let’s first Load all the required libraries

#setwd("C:/GDrive/chingtea/BUDT758T2020/code/timeSeries/") #Don't forget to set your working directory before you start!

library("tidyverse")
library("fpp3")

── Attaching packages ────────────────────────────────────────────────────────────────────────────── fpp3 0.2 ──
✔ tsibble     0.8.6     ✔ feasts      0.1.3
✔ tsibbledata 0.1.0     ✔ fable       0.1.2
── Conflicts ───────────────────────────────────────────────────────────────────────────────── fpp3_conflicts ──
✖ fabletools::accuracy()   masks yardstick::accuracy()
✖ lubridate::date()        masks base::date()
✖ scales::discard()        masks purrr::discard()
✖ plotly::filter()         masks dplyr::filter(), stats::filter()
✖ fabletools::generate()   masks infer::generate()
✖ tsibble::id()            masks dplyr::id()
✖ tsibble::interval()      masks lubridate::interval()
✖ dplyr::lag()             masks stats::lag()
✖ caret::lift()            masks purrr::lift()
✖ fabletools::MAE()        masks caret::MAE()
✖ dials::margin()          masks ggplot2::margin()
✖ tsibble::new_interval()  masks lubridate::new_interval()
✖ fabletools::null_model() masks parsnip::null_model()
✖ fabletools::RMSE()       masks caret::RMSE()

library("plotly")
library("skimr")
library("lubridate")

tsDia <- read_csv('antidiabetic.csv')

Missing column names filled in: 'X1' [1]Parsed with column specification:
cols(
  X1 = col_double(),
  time = col_double(),
  value = col_double()
)

Read the data

month <- as_tibble(yearmonth(seq(as.Date("1991-07-01"), as.Date("2008-06-01"), by = "1 month"))) %>% 
  rename(month=value)


tsDia <- read_csv('antidiabetic.csv') %>%
  bind_cols(month) %>% 
  select(month, drugSales = value) %>%
  as_tsibble(index = month)

Missing column names filled in: 'X1' [1]Parsed with column specification:
cols(
  X1 = col_double(),
  time = col_double(),
  value = col_double()
)

---

Understanding time series data

Boston marathon

Trend without seasonality

tsBoston <-
  boston_marathon %>% 
  filter(Event == "Men's open division") %>%
  select(Year, Time) %>%
  as_tsibble()

Winning times at the Boston Marathon

plotBoston <-
  tsBoston %>% 
  autoplot() +
  xlab("Year (yearly data)") + ylab("Time in HH:MM:SS") +
  ggtitle("Winning times at the Boston Marathon (Men's open division)")

Plot variable not specified, automatically selected `.vars = Time`

ggplotly(plotBoston)

NA

---

Antidiabetic drug sales

Trend and seasonality are both evident

plotDia <-
  tsDia %>% 
  autoplot() +
  xlab("Year") + ylab("Total prescriptions ($ million)") +
  ggtitle("Antidiabetic drug sales") +
  scale_x_date(date_breaks = "years" , date_labels = "%y")
ggplotly(plotDia)

Antidiabetic drug decomposed

Seasonal plot for antidiabetic drug sales

plotDiaSeason <-
  tsDia %>%
  gg_season(drugSales, labels = "both") +
  xlab("Year") + ylab("Total prescriptions ($ million)") +
  ggtitle("Antidiabetic drug sales")
plotDiaSeason

---

U.S. gasoline supplied

Both trend and seasonality

tsGas <-us_gasoline

U.S. finished motor gasoline product supplied

plotGas <-
  tsGas %>% 
  autoplot() +
  xlab("Year (weekly data)") + ylab("Million barrels per week") +
  ggtitle("U.S. finished motor gasoline product supplied") +
  scale_x_date(date_breaks = "years" , date_labels = "%y")
ggplotly(plotGas)

U.S. gasoline supplied decomposed

plotGasDecomposed <- 
  tsGas %>%
  model(STL(Barrels ~ trend() + season(window='periodic'), robust = TRUE)) %>% 
  components() %>%
  autoplot() +
  xlab("Year (weekly data)") + ylab("Million barrels per week") +
  ggtitle("Seasonal and Trend decomposition using Loess (STL) for U.S. gasoline data") +
  scale_x_date(date_breaks = "years" , date_labels = "%y")
ggplotly(plotGasDecomposed)
plotGasDecomposed

---

U.S. retail employment Data

Trend overcomes seasonality

tsEmployment <-
  us_employment %>% 
  filter(Title == "Retail Trade", Month >= '1980-01-01') %>%
  select(Month, Employed)

U.S. retail employment

plotEmp <-
  tsEmployment %>% 
  autoplot() +
  xlab("Year (monthly data)") + ylab("Number of employed in retail (000)") +
  ggtitle("U.S. retail employment data")
ggplotly(plotEmp)

U.S. retail employment decomposed

plotEmpDecomposed <- 
  tsEmployment %>%
  model(STL(Employed ~ trend() + season(window='periodic'), robust = TRUE)) %>% 
  components() %>%
  autoplot() +
  xlab("Year (monthly data)") + ylab("Number of employed in retail (000)") +
  ggtitle("Seasonal and Trend decomposition using Loess (STL) for U.S. retail employment data") +
  scale_x_date(date_breaks = "years" , date_labels = "%y")
plotEmpDecomposed

---

U.S. personal consumption expenditure data

Both trend and seasonality are weak

tsEcon <- us_change %>% 
  filter(Quarter >= '1980-01-01') %>%
  select(Quarter, Consumption)

U.S. personal consumption

plotEcon <-
  tsEcon %>% 
  autoplot() +
  xlab("Year (quarterly data)") + ylab("Percentage changes in personal consumption") +
  ggtitle("U.S. personal consumption expenditure data") +
  scale_x_date(date_breaks = "years" , date_labels = "%y")
ggplotly(plotEcon)

U.S. consumption decomposed

plotEconDecomposed <- 
  tsEcon %>%
  model(STL(Consumption ~ trend() + season(window='periodic'), robust = TRUE)) %>% 
  components() %>%
  autoplot() +
  xlab("Year (quarterly data)") + ylab("Percentage changes in personal consumption") +
  ggtitle("Seasonal and Trend decomposition using Loess (STL) for U.S. consumption data") +
  scale_x_date(date_breaks = "years" , date_labels = "%y")
plotEconDecomposed

---

Accidental deaths in the U.S.

Sesonality without trend

tsAcc <-
  as_tsibble(USAccDeaths) %>% 
  rename(numOfAccDeaths = value)

Number of accidental deaths in the U.S.

plotAcc <-
  tsAcc %>% 
  autoplot() +
  xlab("Year") + ylab("Number of accidental deaths") +
  ggtitle("Number of accidental deaths in the U.S.") +
  scale_x_date(date_breaks = "years" , date_labels = "%y")
ggplotly(plotAcc)

Accidental deaths decomposed

plotAcc <- 
  tsAcc %>%
  model(STL(numOfAccDeaths ~ trend() + season(window='periodic'), robust = TRUE)) %>%
  components() %>%
  autoplot() +
  xlab("Year (monthly data)") + ylab("Number of accidental deaths") +
  ggtitle("Number of accidental deaths in the U.S.") +
  scale_x_date(date_breaks = "years" , date_labels = "%y")
plotAcc

Accidental deaths seasonal subseries

plotAccSub <- 
  tsAcc %>% 
  gg_subseries(numOfAccDeaths) +
  ylab("Number of accidental deaths") +
  xlab("Year (monthly data)") +
  ggtitle("Seasonal subseries plot: Number of accidental deaths in the U.S.")
plotAccSub

---

Preliminary analysis of time series data

Autocorrelation function plot for the antidiabetic drug sales data

plotDiaACF <-
  tsDia %>%
  ACF(drugSales, lag_max = 48) %>%
  autoplot() +  ggtitle("Autocorrelation function (ACF) plot for the antidiabetic drug sales data")
plotDiaACF

Lag plots for the accidental deaths data

plotAccLag <-
  tsAcc %>%
  gg_lag(numOfAccDeaths, geom='point') +
  xlab(NULL) + ylab(NULL) +
  ggtitle("Lag plots for the accidental deaths data")
ggplotly(plotAccLag)

NA

Autocorrelation function (ACF) plot for the accidental deaths data

plotAccACF <-
  tsAcc %>%
  ACF(numOfAccDeaths, lag_max = 9) %>%
  autoplot() +  ggtitle("Autocorrelation function (ACF) plot for the accidental deaths data")
plotAccACF

Noise(randomly generated)

set.seed(333)
y <- tsibble(sample = 1:100, Noise = rnorm(100), index = sample)
y %>%
  autoplot(Noise) + ggtitle("Noise over time") + xlab("Time")

Autocorrection for the noise

y %>%
  ACF(Noise) %>% autoplot()

Additive vs. multiplicative decomposition

Remember antidiabetic drug sales?

plotDia <-
  tsDia %>% 
  autoplot() +
  xlab("Year") + ylab("Total prescriptions ($ million)") +
  ggtitle("Antidiabetic drug sales") +
  scale_x_date(date_breaks = "years" , date_labels = "%y")
ggplotly(plotDia)

Additive decomposition

tsDia %>% 
  model(classical_decomposition(drugSales, type = "additive")) %>%
  components() %>%
  autoplot() +
  xlab("Year") + ylab("Total prescriptions ($ million)") +
  ggtitle("Classical additive decomposition of antidiabetic drug sales") +
  scale_x_date(date_breaks = "years" , date_labels = "%y")

Remember the STL decomposition?

plotDiaDecomposed <- 
  tsDia %>%
  model(STL(drugSales ~ trend(window=10) + season(window='periodic'), robust = TRUE)) %>% 
  components() %>%
  autoplot() +
  xlab("Year") + ylab("Total prescriptions ($ million)") +
  ggtitle("Seasonal and Trend decomposition using Loess (STL decomposition)") +
  scale_x_date(date_breaks = "years" , date_labels = "%y")
plotDiaDecomposed

Multiplicative decomposition

tsDia %>% 
  model(classical_decomposition(drugSales, type = "multiplicative")) %>%
  components() %>%
  autoplot() +
  xlab("Year") + ylab("Total prescriptions ($ million)") +
  ggtitle("Classical multiplicative decomposition of antidiabetic drug sales") +
  scale_x_date(date_breaks = "years" , date_labels = "%y")

---

Decompose data without seasonality(STL)

plotBostonDecomposed <-
  tsBoston %>%
  mutate(Seconds = as.numeric(Time, units="seconds")) %>%
  select(Year, Seconds) %>%
  model(STL(Seconds ~ trend(), robust = TRUE)) %>%
  components() %>%
  autoplot() +
  xlab("Year (yearly data)") + ylab("Time in seconds") +
  ggtitle("Winning times at the Boston Marathon (Men's open division)")
plotBostonDecomposed

Removing seasonal variation(STL)

U.S. retail employment data Seasonally Adjusted

plotEmpSeasonallyAdjusted <-
  tsEmployment %>%
  autoplot(Employed, color='#A9A9B0') +
  autolayer(components(tsEmployment %>% model(STL(Employed))), season_adjust, color='#1490D4') +
  xlab("Year (monthly data)") + ylab("Number of employed in retail (000)") +
  ggtitle("U.S. retail employment data")
ggplotly(plotEmpSeasonallyAdjusted)

us_employment %>%
  filter(Title %in% c('Construction', 'Manufacturing', 'Leisure and Hospitality')) %>% 
  as_tsibble(key=Title) %>%
  select(-Series_ID) %>% 
  features(Employed, feat_stl) %>% 
  rename_all(~str_remove(.x, '_year'))